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Abstract 

The strong-property-fluctuation theory (SPFT) provides a general framework for estimating 
the constitutive parameters of a homogenized composite material (HCM). We developed the 
elastodynamic SPFT for orthotropic HCMs, in order to undertake numerical studies. A specific 
choice of two-point covariance function — which characterizes the distributional statistics of 
the generally ellipsoidal particles that constitute the component materials — was implemented. 
Representative numerical examples revealed that the lowest-order SPFT estimate of the HCM 
stiffness tensor is qualitatively similar to the estimate provided by the Mori-Tanaka mean-field 
formalism, but the differences between the two estimates vary as the orthotropic nature of the 
HCM is accentuated. The second-order SPFT provides a correction to the lowest-order estimate 
of the HCM stiffness tensor and density. The correction, indicating effective dissipation due to 
scattering loss, increases as the HCM becomes less orthotropic but decreases as the correlation 
length becomes smaller. 

Keywords: Homogenization; strong-property-fluctuation theory; metamaterials; Mori-Tanaka 
mean-field formalism. 

1 Introduction 

How do we estimate the effective constitutive properties of composite materials? This question has 
long been considered in the context of acoustics, elastodynamics and electromagnetics [1, 2]. An 
upsurge in interest in this topic has been prompted by the recent proliferation of metamaterials, 
both as theoretical concepts and as physical entities. An operational definition of a metamaterial is 
as an artificial composite material which exhibits properties that are not exhibited by its component 
materials or at least not exhibited to the same extent by its component materials [3]. Metamaterials 
are often exemplified by homogenized composite materials (HCMs). Typically, metamaterials are 
associated with constitutive parameter regimes which have not been accessible conventionally. For 
example, in relation to elastodynamics, metamaterials with negative mass density [4] and negative 
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stiffness [5, 6] have recently been described, whereas negatively-refracting metamaterials have been 
the subject of intense research activity lately in electromagnetics [7]. 

We focus here on the effective elastodynamic properties of a composite material. The HCM 
considered arises in the long-wavelength regime from component materials which are generally 
orthotropic, viscoelastic and randomly distributed as oriented ellipsoidal particles. Our study is 
based on the strong-property-fluctuation theory (SPFT) which — by allowing for higher-order 
characterizations of the distributional statistics of the component materials — provides a multi- 
scattering approach to homogenization [8]. This distinguishes the SPFT from certain well-known 
self consistent approaches to homogenization [9, 10, 11, 12], although we note that more sophis- 
ticated sclf-consistcnt theories have been proposed in recent years [13, 14, 15]. While the general 
character of the SPFT approach to homogenization is reminiscent of multi-scattering theories 
[16, 17, 18], the SPFT provides an estimate of the HCM's constitutive parameters whereas multi- 
scattering approaches generally provide effective wavenumbers [19, 20, 21]. A distinctive feature 
of the SPFT is that it incorporates a renormalized formulation which can accommodate relatively 
strong variations in the constitutive parameters of the component materials. This is because the 
perturbative scheme for averaging the renormalized equations in the SPFT is based on parame- 
ters which remain small even when there are strong fluctuations in the constitutive parameters 
describing the component materials. In contrast, conventional variational methods of homogeniza- 
tion [22, 23, 24, 25, 26] yield bounds which are widely separated when there are large differences 
between the constitutive parameters of the component materials. 

The SPFT has been widely utilized to estimate the electromagnetic constitutive parameters of 
HCMs [27, 28, 29, 30, 31]. Acoustic [32] and elastodynamic [33] versions of the theory have also been 
developed. The general framework for the elastodynamic SPFT, applicable to linear anisotropic 
HCMs, was established in 1999 [33], but no numerical studies have been reported hitherto. In the 
following we apply this theory to examine numerically the case wherein the component materials 
are generally orthotropic materials which are distributed as oriented ellipsoidal particles. Prior to 
undertaking our numerical study, we derive new theoretical results in two areas: 

(i) in the implementation of a two-point covariance function which characterizes the distributions 
of the component materials, and 

(ii) in the simplification of certain integrals in order to make them amenable to numerical com- 
putation. 

The SPFT estimates of the HCM constitutive parameters are illustrated by means of numerical 
examples, and results are compared to those provided by the Mori-Tanaka mean-field approach 
[34, 35]. 

2 Theory 

2.1 Preliminaries 

In applying the elastodynamic SPFT formalism, it is expedient to adopt both matrix and tensor rep- 
resentations [36]. The correspondence between the two representations is described in Appendix A. 
Matrixes are denoted by double underlining and bold font, while vectors are in bold font with no 
underlining. Tensors are represented in normal font with their components indicated by subscripts 
(for nth-order tensors, with < 4) or subscripts and superscripts (for eighth-order tensors). All 
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tensor indexes range from 1 to 3. The pqth component of a matrix A is written as [A] , while 
the pth component of a vector b is written as [h]^. A repeated index implies summation. Thus, we 
have the matrix component [A-B]_^. = [A] [B] vector component [A-b] = [A] [b] , 
and scalar a • b = [aj^ [bj^. The adjoint, determinant and trace of a matrix A are denoted by 
adj ( A ) , det ( A ) and tr ( A ) , respectively. The prefixes Re and Im are used to signify real and 
imaginary parts, respectively, while i = ^/^. 

The SPFT is developed in the frequency domain wherein the stress, strain, and displacement 
have an implicit exp {—iuoi) dependency on time t, oj being the angular frequency. Thus, these are 
generally complex-valued quantities. In order to retrieve the corresponding time-domain quantities, 
the inverse temporal Fourier transform operation must be performed, although one must bear in 
mind that homogenization is essentially a long-wavelength procedure [37, 38]. The possibility of 
viscoelastic behaviour is accommodated through complex-valued constitutive parameters. Stiffness 
tensors are taken to exhibit the usual symmetries 

C'lmpq ~ ^mlpq ~ ^Imqp ~ ^pqlmi (1) 

whilst noting that the symmetry IvaCimpq = ImC'pg/^ has not been proved generally [39]. On 
account of the symmetries (1), the matrix counterpart of tensor Cimpq — namely, the 9x9 stiffness 
matrix C — is symmetric.^ 



2.2 Component materials 

We consider the homogenization of a two-component composite material. The component ma- 
terials, which are themselves homogeneous, are randomly distributed throughout the mixture as 
identically-oriented, conformal, ellipsoidal particles. For convenience, the principal axes of the 
ellipsoidal particles are taken to be aligned with the Cartesian axes. Thus, the surface of each 
ellipsoidal particle may be parameterized by the vector 

r(«) = r?n • r, (2) 
where 77 is a linear measure of size, f is the radial unit vector and the diagonal shape matrix 

U = -7^l 6 1, (a,6,ceR+). (3) 




Let the space occupied by the composite material be denoted by V. It is partitioned into 
parts V^^^ and V^'^^ containing the two component materials labelled as '1' and '2', respectively. 
The distributional statistics of the component materials are described in terms of moments of the 
characteristic functions 

$W(r) = <| (^ = 1,2). (4) 



''Alternatively, in light of (1), the stiflFness tensor may be represented by a symmetric 6x6 matrix [40], but the 
following presentation of the SPFT is more straightforwardly presented in terms of the 9x9 matrix representation. 
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The volume fraction of component material ^, namely /^^^ , is given by the first statistical moment 
of ; i.e., 

($W(r))=/W, (^ = 1,2), (5) 



where the angular brackets denote the ensemble average of the quantity enclosed. Notice that 
/^^^ + /^^^ = 1- The second statistical moment of constitutes a two-point covariance function. 
The physically-motivated form [41] 



($W(r)$W(r')) 



' ($W(r))($W(r'): 
. ($W(r)), 



• (r-r')l > L, 
I • (r - r') I < -L , 



(6) 



is adopted, where L > is the correlation length which is taken to be much smaller than the 

elastodynamic wavelengths but larger than the sizes of the component particles. In the context of 
the electromagnetic SPFT, the specific form of the covariance function has only a secondary influ- 
ence on estimates of HCM constitutive parameters, for a range of physically-plausible covariance 
functions [42]. 

The elastodynamic properties of the component materials '1' and '2' are characterized by their 



stiffness tensors c\}^ and CP^ (or, equivalently, their 9x9 stiffness matrixes C^^\ I G {1, 2}), and 



.(2) 



Impq 



Irnpq 



their densities p^^^ and The stiffness tensors exhibit the symmetries represented in (1). The 
component materials are generally orthotropic [40] in the following developments; i.e., the stiffness 
matrix for each component material may be expressed as 



/ M^^^ 



(^=1,2), 



(7) 



where MS^^ and V^^^ are symmetric and diagonal 3x3 matrixes, respectively, and is the 3x3 null 
matrix. For the degenerate case in which the component material is isotropic, we have 



(8) 





11 


c(^)" 


22 




= aW+2//W 1 

33 






12 




13 




23 


(£=1,2), 




44 




55 




66 > 





where A*^^) and /x*^^) are the Lame constants [43] . 



2.3 Comparison material 

A central concept in the SPFT is that of a homogeneous comparison material. This provides the 
initial ansatz for an iterative procedure that delivers a succession of SPFT estimates of the consti- 
tutive properties of the HCM. As such, the comparison material represents the lowest-order SPFT 
estimate of the HCM. Since we have taken the component materials to be generally orthotropic 
and distributed as ellipsoidal particles, the comparison material is generally orthotropic"'. While 

^In fact, the comparison material would also be orthotropic if (i) the components materials were isotropic but 
distributed as aligned ellipsoidal particles; or (ii) the components materials were orthotropic but distributed as 
spherical particles 
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this is a physically-reasonable assumption here, we remark that the form of the HCM stiffness 
tensor may be derived via certain asymptotic approaches to homogenization [44]. The orthotropic 
comparison material (OCM) is characterized by its stiffness tensor C^^!^ and density with 



^impq exhibiting the symmetries (1). 



The SPFT formulation exploits the spectral Green function of the OCM, which may be expressed 
in 3x3 matrix form as 

g(ocm) (^1^) ^ \eg{k) - io'^p^""'^^-^ ~\ (9) 

with I being the 3x3 identity matrix and a(k) the 3x3 matrix with entries 



[i(k): 



mp 



fe2 



(10) 



Herein, k = fek = (/^i, A;2, fcs) with k = (sin cos ijii, sin^sin^, cos^). For use later on in §2.4, we 
remark that G(°'='") (k) may be conveniently expressed as 



1 



A(k) 



N(k), 



(11) 
(12) 
) • (13) 

A key step in the SPFT — one which facilitates the calculation of c\"^"^ and p^"^"*) — is the 



with the 3x3 matrix function 

^(k) = fe^adj [^(k)] {i(k) - tr [^(k)] l} + (t^V"^"*))^! 



and the scalar function 
A(k) = it^det 



^(k)] -a;^^'"'"^ kh^ {adj [^(k) ] } + (wV^""""^) ^^tr [^(k)' 

Impq 



imposition of the conditions [33, cqs. (2. 72), (2. 73)] 



0, 



^$(i)(r) - + $(2)(r) (p(2) _ p{ocm)y^ ^ 0^ 

in order to remove certain secular terms. In (14), the quantities 

_ ( r<(^) _ Mocm)\ (p _ ^ <)\ 

'ilmpq ~ y^lmst '^Imst ) 'Istpqi l*^ — J-)^J) 



where rfstpq '^^ given implicitly via 

^pq 



IpqstJ st ' 



Jij ~ ^ij ^ '-'vim y-^impq '-'Impq ) ^pq i 



(14) 
(15) 

(16) 

(17) 
(18) 



and the renormalization tensor 



^rstu 



1 

8^ 



2n 



'smy X 



iW' ■ k)t {(n-^ • k).[i-nM-' • k)],„ + (n-^ • \^)r[^-\W' ■ ^)]su} 



(U-1 • k) • (U-i • k) 



(19) 
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Upon substituting (16)-(18) into (14), exploiting (5), and after some algebraic manipulations, 
we obtain 



/ 



(1) 



t 



(20) 



wherein the 9x9 matrix equivalents of the tensors C^^^ and Srstu (namely, C*-°'^™^ and S) have 



been introduced and ^ denotes the matrix operation defined in Appendix A. The OCM stiffness 
matrix may be extracted from (20) as 



g(ocm) ^ g(l) ^ ^{2) (g(2) j^g(l) _ g(2) ^ ^ 



(21) 



where r is the 9x9 matrix representation of the identity tensor Trstui as described in Appendix A. 
This nonlinear relation (21) can be readily solved for C^"*^"*^ by numerical procedures, such as the 
Jacobi method [45]. 

By combining (5) with (15), it follows immediately that the OCM density is the volume average 
of the densities of the component materials '1' and '2'; i.e.. 



piocm) ^ fi^)pW + /(2)p(2). 



(22) 



2.4 Second-order SPFT 

The expressions for the second-order® estimates of the HCM stiffness and density tensors, as derived 
elsewhere [33, eqs. (2.77), (2.78)], are 



Impq Impq o 



d'k^B\^;^{\,) [£(--)(k)]^Jfc, [^-i(k) 



and 



mp 



a-^(k)] } 
(23) 

(24) 



respectively, wherein 5mp is the Kronecker delta function. The eighth-order tensor B|^g*(k) and 
scalar -B(k) represent the spectral covariance functions given as 



5l:r;/(k) 

S(k) 



(2) _ ^ M2) _ .(1) 

Imrs ^Irars J \^tupq >>tupq 



87r3 



j (fR r(R) exp {-ik ■ R) 



('o(2) - f 

^ / (fR r(R) exp (-zk • R) 



(25) 



with 



r(r - r') = ( (r) (r') ) - ( (r) ) ( (r') ) = ( $(2) (r) $(2) (r') ) - ( $(2) (r) ) ( $(2) (r') ) . 

(26) 



The first-order SPFT estimate is identical to the zeroth-order SPFT estimate which is represented by the 
comparison material. 
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We now proceed to simplify the expressions for c\^^^ and p^^^^ presented in (23) and (24), 
in order to make them numericahy tractable. Wc begin with the integral on the right sides of (25) 
which, upon implementing the step function-shaped covariance function (6), may be expressed as 



I SR r(R) exp (-zk • R) = /" d^R exp [-i (U • k) • R] 

J J\R.\<L ~ 

Thus, we find that Bl'!^pg(k) and -B(k) are given by 

/■(I) f (2) /'t(2) A /^t(2) ^ r . / 

•I J y'lmrs ^ImrsJ yitupq Stupg J gin (kaL) 



(27) 



5!«p7(k) 
B(k) 



/(l)/(2) (p(2)_^(l)) 



2 {-Kka) 

(1)n2 



L cos {kaL) 



2 (tt/cct)^ 
wherein the scalar function 

a = a{e,(t)] 



sin (kaL) 
ka 



Lcos (kaL) 



(28) 



sin^ cos^ d) + b"^ sin^ ^ sin^ + cos^ 0. 



(29) 



Upon substituting (28) into (23) and (24), the integrals therein with respect to k can be evalu- 
ated by means of calculus of residues: The roots of A(k) = give rise to six poles in the complex-fc 
plane, located at k = ±pi, ±p2 and ±^3, chosen such that Re pi > {i = 1,2,3). From (13), we 
find that. 



Pi = Pa--^ 



2V3p, 



B 



Pc 



Pc det 



a(k) 



2V3det 



a(k) 



pI 



Pa + 



Pa + 



;i + /V3)p« 



(i-/:v^)Pc 



22/3p^det ^(k) 

(1 - iV3)PB 



24/3 det 



a(k) 



(1 + iV^)Pc 



22/3Pcdet a(k) 



24/3 det 



a(k) 



(30) 
(31) 
(32) 



wherein 



Pa 
Pb 

p3 





adj 


i(k)]} 


3 det 


i(k)^ 





2 / ' 

^^2^{ocm)^ (3 det [^(k)] tr [^(k)] -trjadj [i(k)]}^ 



= Po + ^^P^ + Ph, 
Pd = (a;V°^"^))^ (2tr{adj [^(k)]}' 



-9 det 



^(k)] trjadj ^(k)] } tr [^(k) +27 det ^(k) ^) 



(33) 

(34) 
(35) 

(36) 
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Thus, by application of the Cauchy residue theorem [46], the SPFT estimates are dehvered as 



(spft) 
Irnpq 



. ,2 Jocm) f{l) f(2) ( pm Wt(2) ^ 

^'Impq ^ '■ '■ ^ 



27r 



A-ni 



d(j)de- 



J 0=0 



kt sin 9 




rv 






{kaf det 







(37) 



and 



^2^(1)^(2) (^(2)_^(l)y 

27ri 



I / ' 

J 6=0 J 9=0 



>de 



smt 



det 



a(k) 



b(k) 



mp 



(38) 



where 



b(k) = 



1 r e^-^'^Pi^(pin • k) 
giLap3N(p3U • k) 



1 — iLapi 



1 — iLaps 



giLap2N(p2U • k) 



1(0) ' 



pi) (pi - pi) 



1 — iLap2 



(39) 



(^pi(pi-pi)(pf-pi) 

The integrals in (37) and (38) are readily evaluated by standard numerical methods [47]. 

Significantly, the second-order SPFT estimates C^^^^^ and p^mp^*^ are complex-valued even 

when the corresponding quantities for the component materials, i.e., C^^^^ and p^^-* {^ = 1,2), 
are real-valued. This reflects the fact that the SPFT effectively takes into account losses due to 
scattering. This feature is not unique to the SPFT: it arises generally in multi-scattering approaches 
to homogenization [16, 20, 48]. We note that for 



{spft) 



is required to be 



(i) the timc-avcragcd strain energy density to be positive-valued, ReC 
positive-definite; and 

(ii) the time-averaged dissipated energy density to be positive-valued, — ImC is required to 
be positive-semi-definite, 



" (spft) 

where C is the 6x6 matrix with components 



^{spft) 



St 



Q{spft) 



St 



{s,t G {1,2,. ..,6}) 



and ^^Pf^^ is the 9x9 matrix equivalent to the SPFT stiffness tensor Ci^^^ . 

It is notable too that the second-order SPFT estimates c\^^^ and Pmp'^*'' are explicitly depen- 
dent on frequency, whereas the corresponding zeroth-order SPFT estimates exhibit only an implicit 
dependency on frequency via the frequency-dependent constitutive parameters of the component 
materials. Accordingly, the second-order SPFT estimates may be viewed as low frequency correc- 
tions to the quasi-static estimates provided by the the zeroth-order SPFT. 

A complex-valued anisotropic density, as delivered by (38), is not without precedent [49]; see 
Milton [50] for a discussion on this issue. 
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3 Numerical results 



Let us now illustrate the theory presented in §2 by means of some representative numerical examples. 
We consider homogenizations wherein the component materials are acetal and glass (or orthotropic 
perturbations of these in §3.2). The corresponding results are qualitatively similar to those we 
found from homogenizations involving a wide range of different component materials, characterized 
by widely different constitutive parameters, which are not presented here. In order to provide a 
baseline for the SPFT estimate of the HCM stiffness tensor, the corresponding results provided by 
the Mori-Tanaka mean-field formalism [34, 35] were also computed. The Mori-Tanaka formalism 
was chosen as a comparison for the SPFT because it is well-established and straightforwardly 
implemented [36, 51]. Comparative studies involving the Mori-Tanaka and other homogenization 
formalisms are reported elsewhere; see [52, 53, 54], for example. The Mori-Tanaka estimate of the 
9x9 stiffness matrix of the HCM may be written as [51] 



t 



(40) 



where 



(41) 



and s'^^*'^^ is the 9x9 Eshelby matrix [55]. The evaluation of this matrix is described in Appendix B. 

In the remainder of this section, we present the 9x9 stiffness matrix of the HCM, namely 
Q{hcm)^ as estimated by the lowest-order SPFT (i.e., hem = ocm), the second-order SPFT (i.e., 
hem = spft) and the Mori-Tanaka mean-field formalism (i.e., hem = MT). The matrix Qi^^"^) 
generally has the orthotropic form represented in (7) with i = hem. We also present the second- 
order SPFT density tensor pmp^^^; numerical results for the lowest-order SPFT density p^""^"*) 
need not be presented here as that quantity is simply the volume average of the densities of the 
component materials. For all second-order SPFT computations, we selected a; = 27r x 10^ s~^. 



3.1 Isotropic component materials distributed as oriented ellipsoidal peirticles 

Let us begin by considering the scenario in which the component materials are both isotropic. The 
component material '1' is taken to be acetal (i.e., A^^^ = A^'*'^'^^ fi^^^ = /i'^"'^'^) and p^^^ = p'-"^'^^), and 
component material '2' to be glass (i.e., A*^^) = A*^^'"^, /x^^^ = /x^^'") and p^^^ = p^^''")). The Lame 
constants and densities for these two materials are as follows [56, 57]: 



^(ace) 



2.68 GPa, = 1.15 GPa, 



21.73 GPa, fi 



(gla) 



29.2 GPa, p 



piace) = 1,43 X kgm-''^ 

(gla) ^ 2.23 X 10^ kgm-3 



(42) 



The eccentricities of the ellipsoidal component particles are specified by the parameters {a,b,e}, 

per (2) and (3). 

In Fig. 1 the components of the HCM stiffness matrix q(^^"^) ^ as computed using the lowest- 
order SPFT and the Mori-Tanaka formalism, are plotted as functions of volume fraction /^^^ for 



the case a = b = e. Since the HCM is isotropic in this case, only the components 



Qf{hcTn) 



11 



12 



44 



^(hcm) 

are presented, per (8) with 



9 



i = hem. Notice that the following limits necessarily apply for both the SPFT and Mori-Tanaka 
estimates: 

lim C('^"™) = CW, lim C^hcm)^Q{2)_ U2,) 

/(2)_0= = /(2)-^l = 

I : is apparent from Fig. 1 that, while the lowest-order SPFT and the Mori-Tanaka estimates are 
qualitatively similar, the Mori-Tanaka estimates display a greater deviation from the naive HCM 
estimate 

/(I) c(^) + C(2) for mid-range values of /^^^. For further comparison in this 



pq 



pq 



isotropic scenario, the familiar variational bounds on QK^'^i^) and QV^^-"'-) established by 

L= Jii L= J44 

Hashin and Shtrikman [2, 22] are also presented in Fig. 1: the lower Hashin Shtrikman bound 
coincides with the Mori-Tanaka estimate and the lowest-order SPFT estimate lies within the 
upper and lower Hashin-Shtrikman bounds for all values of f^'^\ Parenthetically, we note that for 
isotropic HCMs the lowest order SPFT estimates are the same as those provided by the well-known 
formalisms of Hill [9] and Budiansky [11], as demonstrated elsewhere [33]. 

The corresponding lowest-order SPFT and Mori-Tanaka estimates for the orthotropic HCM 
arising from the distribution of component material as ellipsoids described by {a/c = 5, b/c = 1.5} 

are presented in Fig. 2. The matrix entries C'^'*"^'"^ are plotted against /^^^ for pq G {11, 12, 44}. 

L Jpg 

The graphs for pq = 13 and 23 are qualitatively similar to those for pq = 12; and those for pq = 55 
and 66 are qualitatively similar to those for pq = 44. The degree of orthotropy exhibited by the 



HCM can be gauged by relative differences in the values of 



Q(/icm) 



similarly by relative differences in 



Q {hem) 



pq 



for pq e {11, 22, 33} (and 



pq 



for pq G {44, 55, 66} and by relative differences in 



pq 



for pq e {12, 13, 23}). These relative differences are greatest for mid-range values of 



the volume fraction f^'^\ 

The orthotropic nature of the HCM is accentuated by using component materials with more 
eccentrically-shaped particles. This is illustrated by Fig. 3, which shows results computed for 
the same scenario as for Fig. 2 but with ellipsoidal particles described by {a/c = 10, b/c = 2}. A 
comparison of Figs. 1-3 reveals that differences between the estimates of the lowest-order SPFT 
and the Mori-Tanaka mean-field formalism vary slightly as the orthotropic nature of the HCM is 
accentuated. For example, the difference between the lowest-order SPFT and the Mori-Tanaka 



estimates of the 



Qi{hcm) 



44 



increases as the HCM becomes more orthotropic. 



Now let us turn to the second-order SPFT estimates of the HCM constitutive parameters. We 
considered these quantities as functions of kL, where 




,(1) 



+ 



(1) 



,(2) 



+ 



9(2) 



(2) 



(44) 



is an approximate wavenumber calculated as the average of the shear and longitudinal wavenumbers 
in the component materials, and L is the correlation length associated with the the two-point 
covariance function (6). Since L is required to be smaller than characteristic wavelengths in the 
HCM (but larger than the sizes of the component particles), we restrict our attention to < kL < 

0.6. Fig. 4 shows the real and imaginary parts of the components of C * = C^^Pf^^f — 
plotted against kL for /^^^^ = 0.5. The values of the shape parameters {o, b, c} correspond to 
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are 



those used in the calculations for Figs. 1-3. As previously, only the matrix entries 

presented for pq G {11, 12, 44}. The graphs for pq = 13 and 23 are qualitatively similar to those for 
pq = 12; and those for pq = 55 and 66 are qualitatively similar to those for pq = 44. Notice that 



lim C^'Pf^^ 

L^0 = 



and 









< 












pq 








pq 



(45) 



(46) 



for all nonzero matrix entries. Furthermore, for the particular example considered here, the mag- 



nitude of 



C 



(spft) 



pq 



generally decreases as the particles of the component materials become more 



eccentric in shape. 

A very striking feature of the second-order SPFT estimates presented in Fig. 4 is that 



Im 



Qispft) 



pq 



(47) 



whereas Im 


Q{a,b) 


= Im 


Qiocm) 






pq 





pq 



= 0. Furthermore, the magnitude of Im 



Qispft) 



pq 



grows steadily as the correlation length increases from zero. These observations may be interpreted 
in terms of effective losses due to scattering as follows. For all reported calculations, ReC is 

positive definite and — Im C*" ^'^ ^ is positivc-scmi -definite, which together imply that the associated 
time -averaged strain energy and dissipated energy densities arc positive - valued [39], as discussed in 

§2.4. Accordingly, the emergence of nonzero imaginary parts of C(*^-/'*) indicates that the HCM 

L ipq 

has acquired an effectively dissipative nature, despite the component materials being nondissipative. 
This effective dissipation must be a scattering loss, because the second-order SPFT accommodates 
interactions between spatially-distinct scattering centres via the two-point covariance function (6). 
As the correlation length increases, the number of scattering centres that can mutually interact 
also increases, thereby increasing the scattering loss per unit volume. 

Lastly in this subsection, the real and imaginary parts of the second-order SPFT density tensor 
Ppq^'^^ = Ppq'^^^ — p^°'^"^^ are plotted as functions of kL in Fig. 5. Only the p = q components are 
presented, as the p ^ q components are negligibly small. The density tensor exhibits characteristics 
similar to those of the corresponding stiffness tensor insofar as 



lim o^**'-^*) 



and 



zispft) 
rpq 



(spft) 



(ocm) 



(ocm) 



(48) 



(49) 



for all values of the indexes p and q. Also, \Ppq '\ generally decreases as the shape of the particles 
of the component materials deviates further from spherical. 
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3.2 Orthotropic component materials distributed as spheres 

Let us now explore the scenario wherein the component materials are orthotropic perturbations of 
the isotropic component materials considered in §3.1. In the notation of (7), we choose 



and 



(^^{ace) ^ 2/x("'^^)) (1 + q) X^""^^ (1 - <;) 



A(»^e) (1 + 2(^) 



Xiace) (1 + 2^) 
(A(«ce) + 2/x("'=^)) (1 
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(^(ace)^ (1 


















(ace) 



(1 - 1^) 



(50) 



p(2) 



(A(s'") + (1 + 2^ 

A(9'«) (1 - 2q) 

(^W) (1 - §0 



A(ff'") (1-2?) 

(aW + 2/v,(9''^)) (1 + I?) 

A(s'«) (1 - i,j) 



;i + i.) 



(A(9^«) + 2/Lt(^^^")) (l 



(gla) 



(1 - 1^) 








(gla) 



(1 - 10 



(51) 

where the real-valued scalar ? controls the degree of orthotropy. Thus, at fixed values of <; the 

component materials may be viewed as being locally orthotropic. As in §3.1, the densities of the 



.(2) 



(gla) 



The component materials are 



component materials are taken to be p^^^ = p^"^'^) and 
distributed as spherical particles (i.e., a = b = c). 

The lowest-order SPFT and Mori-Tanaka estimates for the HCM arising from orthotropic com- 
ponent materials characterized by ? = 0.05 and ? = 0.1 are presented in Fig. 6 and 7, respectively. 
As previously, only a representative selection of the entries of Q^^^"^) are provided here. The plots 
for = 0, for which case the HCM is isotropic, are the ones displayed in Fig. 1. In a manner 
resembling the scenario considered in §3.1, the lowest-order SPFT and the Mori-Tanaka estimates 
are qualitatively similar, but the Mori-Tanaka estimates display a greater deviation from the naive 



HCM estimate /(^^ 



c(i) 



+ /(2) 



pq 



C(2) 



for mid range values of f^'^\ at all values of ■ 



pq 



The degree of orthotropy exhibited by the HCM clearly increases as the value of increases, 
and differences between the estimates of the lowest-order SPFT and the Mori-Tanaka mean- 
field formalism also vary as ? increases. To explore this matter further, in Fig. 8 the associated 



ratios 



Qi{hcm) 



11 



/ 



Q(/icm) 



33 



12 



/ 



Q {hem) 



and 



23 



{hem) 



44 



/ 



Q {hem) 



44 



are plot- 



ted against /^^^ for ? = 0.05 and 0.1. The three different patterns arc portrayed in the three 



plots: for 



Qf^{hcm) 



11 



/ 



Qf{hcm) 



differences between the lowest-order SPFT and the Mori- 



33 



Tanaka estimates are larger for when the HCM is more orthotropic; the reverse is the case for 



^{hcm) 



12 



I 



Qijyhcrn) 



while for 



23 



Qi{hcm) 



44 



/ 



^{hcm) 



66 



there is no noticeable difference be- 



tween the lowest-order SPFT and Mori-Tanaka estimates as the degree of HCM orthotropy is 
increased. 

Next we focus on the second-order SPFT estimate of the HCM stiffness tensor. The real and 
imaginary parts of a representative selection of entries of ^^'^^^^ = Qi^pf^) _ Qi"cm) graphed 
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against kL in Fig. 9. The volume fraction is fixed at f^"^^ = 0.5. The values of the orthotropy 
parameter ^ are 0, 0.05 and 0.1, in correspondence with the calculations of Figs. 1, 6 and 7. As 

- (spft) 

we observed in §3.1, the magnitude of the components of C generally decrease as the HCM 
becomes more orthotropic. Also, the second-order SPFT estimate C^^Pf^") has components with 
nonzero imaginary parts, which implies that the HCM is effectively dissipative even though the 
component materials are nondissipative. Furthermore, the HCM becomes increasingly dissipative 
as the correlation length increases, this effective dissipation being attributable to scattering losses. 

Finally, the real and imaginary parts of the second-order SPFT density tensor Ppq^^^ = Ppq^^^ — 
p{ocm) g^j.g plotted as functions of kL in Fig. 10. As previously in §3.1, the components ior p ^ q are 
negligibly small so only the p = q components are provided here. The density plots resemble those 
of the corresponding stiffness tensor; i.e., the components p^^*^ are much smaller than p^°'^^) and 
they increase rapidly from zero as L increases from zero. The magnitudes of p^p^^^ are smallest 
when the orthotropy parameter describing the component materials is greatest. 

4 Closing remarks 

The elastodynamic SPFT has been further developed, in order to undertake numerical studies 
based on a specific choice of two-point covariance function. From our theoretical considerations 
in §2 and our representative numerical studies in §3, involving generally orthotropic component 
materials which are distributed as oriented ellipsoids, the following conclusions were drawn: 

• The lowest-order SPFT estimate of the HCM stiffness tensor is qualitatively similar to that 
provided by the Mori— Tanaka mean— field formalism. 

• Differences between the estimates of the lowest-order SPFT and the Mori-Tanaka mean-field 
formalism are greatest for mid-range values of the volume fraction. 

• Differences between the estimates of the lowest-order SPFT and the Mori-Tanaka mean-field 
formalism vary as the HCM becomes more orthotropic. The degree of orthotropy of the HCM 
may be increased by increasing either the degree of orthotropy of component materials or the 
degree of eccentricity (nonsphericity) of the component particles. 

• The second-order SPFT provides a low-frequency correction to the quasi-static lowest-order 
estimates of the HCM stiffness tensor and density. The correction vanishes as the correlation 
length tends to zero. 

• The correction provided by second-order SPFT, though relatively small in magnitude, is 
highly significant as it indicates effective dissipation due to scattering loss. 

• Differences between the second-order and lowest-order SPFT estimates of the HCM consti- 
tutive parameters diminish as the HCM becomes more orthotropic. 

The ability to accommodate higher-order descriptions of the distributional statistics of the 
component materials bodes well for the future deployment of the SPFT in exploring the complex 
behaviour of metamaterials as HCMs. Additionally, since the SPFT has been now established for 

acoustic, electromagnetic and elastodynamic homogenization scenarios, the prospect of considering 
HCMs with mixed acoustic/elastodynamic/electromagnetic constitutive relations beckons. 
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Appendix A 
Matrix/tensor algebra 

A fourth-order tensor A^stu {r,s,t,u G {1,2,3}) has 81 components. If it obeys the symmetries 
A.fstu = Asrtu = Arsut = Aturs, it Can be represented by a 9x9 matrix A with components [ A]^^ 
{R, S G {1,...,9}) . Similarly, the nine entries of a second-order tensor 5,.^ (r, ,s G {1,2,3}) may 
be expressed as a column 9-vector B with components [B]^ (R G {1, . . . ,9}). The scheme for 
converting between the tensor subscript pairs rs and tu and the matrix indexes RS or vector index 
R is provided in Table 1. 



R,S 


rs, tu 


R,S 


rs, tu 


R,S 


rs, tu 


1 


11 


4 


23 or 32 


7 


23 or 32 


2 


22 


5 


13 or 31 


8 


13 or 31 


3 


33 


6 


12 or 21 


9 


12 or 21 



Table 1: Conversion between tensor and matrix/ vector subscripts. 

The most general 9x9 matrix A considered in this paper has the form 




(52) 



where a is a general 3x3 matrix, /3 is a diagonal 3x3 matrix, and is the null 3x3 matrix. If we 
define a 9x9 matrix A''' as [36] 



At 







g 

then A^ • A = A • A^ = r, where r is the 9x9 matrix counterpart of the identity tensor 



1 «-i 1 

4 : 



(53) 



Trstu — 2 i^rt^su "I" ^ru^st) 



(54) 
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Appendix B 

Eshelby matrix/tensor 

If the component materials are orthotropic and distributed as spherical particles (i.e., a = b = c) 
then the tensor counterpart of the 9x9 Eshelby matrix is given as [58] 



1 /""hi /'27r 



(55) 



wherein 



D 



(1) 



D = €mnlKmlKn2Ku, K^, = ClAi&j'&l 



^1 



Ci 



^2 



^3 



(56) 



a 6 ' c 

Ci = (1 - Cl)'/' cos(a;), C2 = (1 - Cl)'/' sin(a;), Cs = Cs J 

with Cijk being the Levi-Civita symbol. The integrals in (55) can be evaluated using standard 
numerical methods [47]. 

If the component materials are isotropic and distributed as ellipsoidal particles described by the 
shape matrix U, then the Eshelby matrix has the form represented in (52) with distinct components 
given as [51] 



g(es/i) 



g(es/i) 



(g(es/i) 



g(es/t) 



g(es/i) 



g(es?i) 



g(es/t) 



11 



12 



13 



21 



22 



23 



31 



Sa^Iaa + (1 


- 2z/ 




87r(l- 


uW] 




SPlafS - (1 - 


- 2u 




87r(l- 






3c laj (l 


- 2i/< 




87r(l- 






Sa'^Iafs - (1 


- 2u 




87r (1 - 


uW] 




362/^^ + (1 - 


- 2v 




87r(l- 






Sc'^Ip^ + (1 - 


- 2i^( 




87r(l- 


uW] 




Sa^Iaj - (1 


- 2u 




87r(l- 







(57) 



(58) 



(59) 



(60) 



(61) 



(62) 



(63) 
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g(es/i) 


36%^ -(1 

QO fl'rr ( 1 


Z/v / 1 






3c2/^^ + (l 


7/(1)^ ' 

i/y^ j 






3(62 ^ ^2^j^^ ^ 


-2z.«)(/^ + /^) 


A A 
J 44 


IDTT I i - 




c* (esfi) 




3(a2 + c2)/„^ + (1 


-2iyW){Ic. + I^) 


55 


167r (l - 


-1.(1)) 


g(es/i) 




3(a2 + 52)/,^ + (1 


-2z/W)(/„ + I^) 


66 


167r (l - 


-1/(1)) 



where i/(i) 
we have 



Ad) 



2(A(i)+/x(i)) 



(64) 



(65) 



(66) 



(67) 



(68) 



is the Poisson ratio of component material '1'. For the case a> b> c 



Anabc 

°= (a2-62)(a2_c2)i/2 
Awabc 



(62 -c2)(a2 -c2)i/2 [ac 



F{e, k) - E{e, k) 
b 



.{a'-cY^-E{e,k) 



Ip = 477 - {la + I-y), Ial3 = 

3^ 



la - 1/3 



let I'y 



I Pi = 



1(3 - I-r 



} , (69) 



Att Att Ait 

Iaa = 17-2- {Ial3 + lafi), 1(3(3 = 177o - i^ap + Ifi-y), -^77 = ^ " ("^"7 + ^/^t) 



3(62 -a2)' 3(c2_o2)' ^/^T 3(c2-62) 

47r 

3^2 ^""P ' "P7)' -'77 ~ 

with the elhptic integrals given by 

re 

E{9, k)= d(t>{l- fc2 sin2 0)i/2 
Jo 
l-e 

F{9, k)= d(t>{l - P sin2 <^)-i/2 
Jo 



(70) 



wherein 
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Figure 1: Plots of 
order SPFT 
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and 
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I.e., 



hem 



ocm] 



(in GPa) as estimated using the lowest- 

(red, solid curves) and the Mori-Tanaka mean-field formalism 
(i.e., hem = MT) (black, dashed curves), against the volume fraction of component material '2'. 
Also plotted are the upper and lower Hashin-Shtrikman bounds (blue, long dashed curves) for 



and 



11 



Qi{hcm) 



the lower Hashin-Shtrikman bounds coincide with the Mori-Tanaka 



44 



estimates. Component material '1' is acetal and component material '2' is glass, as specified in 
(42). The component materials are distributed as spheres (i.e., a = b = c). 
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Figure 2: Plots of 



Qihcm) ^ ^^^Yi rs £ {11,12,44} (in GPa) as estimated using the lowest- 

J rs 

order SPFT (i.e., hem = ocm) (red, solid curves) and the Mori-Tanaka mean-field formalism 
(i.e., hem = MT) (black, dashed curves), against the volume fraction of component medium '2'. 
Component material '1' is acetal and component material '2' is glass, as specified in (42). The 
component materials are distributed as ellipsoids with a/c = 5 and b/c = 1.5. 
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Figure 4: The real and imaginary parts of 



with rs £ {11,12,44} (in GPa), plotted 



as functions of kL, for /^^^ = 0.5. Component medium '1' is acetal and component medium '2' is 
glass, as specified in (42). The component materials are distributed as (i) spheres (i.e., a = b = c) 
(red, solid curves), or (ii) ellipsoids with shape parameters {a/c = 5, b/c = 1.5} (blue, short-dashed 
curves), or (iii) ellipsoids with shape parameters {a/c = 10, b/c = 2} (black, long-dashed curves). 
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Figure 5: As Fig. 4 but the quantities plotted are the real and imaginary parts of the excess of 
be second-order SPFT density tensor oa 
_ p(ocm)^ ^ {1,2,3}), in kgm-3. 



(sv ft) 

the second-order SPFT density tensor over the density of the comparison material, i.e., pr/ 
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Figure 6: Plots of 



Q{hcm) ^ ^-^j^ g {11,12,33,44} (in GPa) as estimated using the lowest- 

J rs 

order SPFT (i.e., hem = ocm) (red, solid curves) and the Mori-Tanaka mean-field formalism (i.e., 
hem = MT) (black, dashed curves), against the volume fraction of component material '2'. The 
component materials are distributed as spheres. Their constitutive parameters are specified by (50) 
and (51), with the orthotropy parameter q = 0.05. 
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Figure 7: As Fig. 6 but with orthotropy parameter = 0.1. 
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Figure 8: Plot of 
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Qi{hcm) j ^(hcm) Q\{hcra) j ^(hcm) and (^(hcm) i Qi{hcm) 

= J 11 L= J33' L= J 12 L= J23 L= J44 L = 

(in GPa) as estimated using the lowest-order SPFT (i.e., hem = ocm) (red, solid curves), the Mori 
Tanaka mean-field formalism (i.e., hem = MT) (black, dashed curves) against the volume fraction 
of component material '2'. Component material '1' is acetal and component material '2' is glass, 
as specified in (42). The component materials are distributed as spheres with the orthotropy 
parameter = 0.05 (thin curves) and <f = 0.1 (thick curves). 
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Figure 9: The real and imaginary parts of 



, with rs £ {11,12,44} (in GPa) plotted 



as functions of kL, for /^^^ = 0.5. The component materials are distributed as spheres. Their 
constitutive parameters are specified by (50) and (51), with the orthotropy parameter ? = (red, 
solid curves), <j = 0.05 (blue, short-dashed curves) and ? = 0.1 (black, long-dashed curves). 
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Figure 10: As Fig. 9 but the quantities plotted are the real and imaginary parts of the excess of 
be second-order SPFT density tensor oa 
_ p(ocm)^ ^ {1,2,3}), in kgm-3. 



the second-order SPFT density tensor over the density of the comparison material, i.e., fxfr^^^ 
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